/*
Replication Code for "When Does Oil Harm Child Mortality?
Authors: Nisha Bellinger and Matthew Fails
Forthcoming in Political Research Quarterly
*/

**************************
**************************

* Read in the data and do some basic data management

use "c:\data\BellingerFailsPRQReplicationData.dta", clear

**Tsset to declare the structure of the data

tsset cowcode year

**Sort the data for variable generation
sort cowcode year

**************************
**************************

*Generate some key variables

*Oil Income and Oil State Status

gen oilincome = ln(oil_gas_valuepop_2014+1)
label var oilincome "log of oil income per capita, newest Ross/Mahdavi data"
sum oilincome
hist oilincome, freq

*generate an oil state dummy variable, where 1 = states who averaged more than $100 per person
bysort cowcode: egen oil_avg = mean(oil_gas_valuepop_2014)
tab oil_avg
bysort cowcode: egen oilinc_avg = mean(oilincome)
bysort cowcode: egen oilinc_8014 = mean(oilincome) if year>=1980 & year<=2014
gen oilstate = 0
replace oilstate = 1 if oil_avg>100 & oil_avg!=.
list CountryName if oilstate==1
label var oilstate "Oil State Dummy, defined as averaging more than $100USD per capita over panel"

gen oilstate_2 = 0 
replace oilstate_2 = 1 if oil_gas_valuepop_2014>100
label var oilstate_2 "Oil State Dummy, define as annual oil income more than $100 USD per capita"

*Generate other control variables
gen lincome_new = ln(WBincome_new)
label var lincome_new "Log income per capita, newest WDI"
*note that this is post 1990 only
*preferred income per capita measure, drawing from Graham IPE with 1999 correction
gen lgdppc = ln(gdppc_wdi)
label var lgdppc "log of gdp per capita, Graham/WDI data"
rename lngdppc_wdi_pw lincome_comb 
gen lincome = ln(income)
label var lincome "Log income per capita, WDI"
gen inc_comb = lincome_comb
replace inc_comb = lincome_new if year==1999
label var inc_comb "Graham WDI PW log income per capita constant, with 2000 replace"
*Graham is missing 1999 value, so replace with my own WDI version of this so I can use the pre 1991 values, 
*but also won't lose 1999 (year 2000 in models with one year lag on IVs)

gen log_pop = ln(population)
label var log_pop "Log population"

**Generate a measure of internal conflict; draws on ucdp type 3 and 4 (dummies for each) from Quality of Government Large Dataset
*this is a measure of internal and internationalized internal conflict, from the UCDP/PRIO datset
gen conflict = 0
replace conflict = 1 if ucdp_type3>0 & ucdp_type3!=.
replace conflict = 1 if ucdp_type4>0 & ucdp_type4!=.
tab conflict
tab ucdp_type3
tab ucdp_type4
label var conflict "Internal Conflict Year dummy variable"

***Generate the dependent variable, based on Child Mortality rates
*DV using natural log of child mortality, to address skewness
gen ln_cm = ln(ChildMortality)
sum ln_cm

**Generate measure of change in CM from 1980-2014 for Figure 1
*note: drop Dominica from sample because of extreme negative value which distorts Y axis
gen cm80 = ln_cm if year==1980
harmby cm80, by(cowcode)
gen cm14 = ln_cm if year==2014
harmby cm14, by(cowcode)
bysort cowcode: gen cm8014 = (cm14-cm80)*-1 if country!="Dominica"
label var cm8014 "Change in natural log of child mortality, 1980-2014"


*************************
**************************

* Create the key motivating graph; our version of Ross' (2012) figure 6.6

*Start by identifying 'long term producers'

*using 'oil80dummy' which is 1 if oilincome>100 for every year from 1980 onward
gen oil100 = 0 if year<=2014
replace oil100 = 1 if oil_gas_valuepop_2014>100 & oil_gas_valuepop_2014!=. & year<=2014
label var oil100 "Dummy variable for whether or not oil income was greater than 100 in that year"

by cowcode: egen oil80 = mean(oil100) if year>=1980
by cowcode: gen oil80dummy = 0
replace oil80dummy = 1 if oil80==1
label var oil80dummy "Dummy variable for all countries which had at least $100 oil income per capita every year since 1980"

*Generate Figure 1 - motivating scatterplot
*Includes three linear fit lines: full sample, long term producer sample, long term producer only autocracy sample	
graph twoway (scatter cm8014 oilinc_8014 if year==2014) ///
	(scatter cm8014 oilinc_8014 if year==2014 & oil80dummy==1, mlabel(country) ///
	scheme(plotplain) legend(off) ylabel(0(1)3) yscale(range(0 3)) ///
	ytitle("Drop in Child Mortality (ln), 1980-2014") ///
	xtitle("Average Oil income per capita (ln), 1980-2014") ///
	xsize(8) note("Labeled countries exceeded $100 oil income per capita for every year 1980-2014")) ///
	(lfit cm8014 oilinc_8014 if year==2014, lpattern(dash)) (lfit cm8014 oilinc_8014 if year==2014 & oil80dummy==1, lpattern(line) lwidth(medthick)) ///
	(lfit cm8014 oilinc_8014 if year==2010 & oil80dummy==1 & gwf_nonautocracy!="democracy", lpattern(line) lwidth(medthin) ///
	note("Labeled countries exceeded $100 oil income per capita for every year 1980-2014. Dashed line is linear fit for full sample." ///
	"Thick solid line is linear fit for labeled country sample. Thin solid line is linear fit for labeled country sample excluding democracies.") ///
	title("Figure 1: Oil Income and Reductions in Child Mortality, 1980-2014" " ", pos(11)) ///
	)
	
*further evidence of variation within non-democracies
sum cm8014 if gwf_nonautocracy!="democracy" & year==2010
sum cm8014 if gwf_nonautocracy=="democracy" & year==2010

	
	
**************************
**************************

* Start Regression Analysis

*Estimate 'standard' model (LDV, unit and time fe's)
		
by cowcode: gen ldv = l.ln_cm
*must generate a lagged dv as a separate variable for the 'margins' procedure, and not use the factor notation

*estimate model and store results
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
estimates store m1

*Plot these regression coefficients as Figure 2
coefplot m1, drop(_cons *.year ldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income") ///
	title("Figure 2: Fixed Effect Regression Coefficients and 95% and 90% CI's" " ", pos(11)) ///
	xsize(8) note(" " "DV is (log) Child Mortality. Sample includes all autocratic regime years 1980-2010. All IV's lagged one year." ///
	"Lagged DV, constant, year, and country fixed effects included but not reported. N = 2,135, Clusters = 103." ///
	"Standard errors clustered on country. See text for full sample, variable, and model description." ///
		)

*Calculate and plot marginal effects
*Figure 3 - marginal effect of oil income on child mortality, across range of autocratic personalization scores
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%", axis(2)) ///
	ytitle("Distribution of autocratic personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(8) ///
	title("Figure 3: Marginal Effect of Oil Income on Child Mortality" " ", pos(11))
	
*Appendix Figure XXXX - symmetrical interpretation of interaction effect	
margins, dydx(l.latent_personalism) cont at(l.oilincome = (0(1)11)) atmeans
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Oil income per capita (log)") ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist oilincome if e(sample), ///
	blcolor(gs10) fcolor(white) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%" , axis(2)) ///
	ytitle("Distribution of oil income, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(8) ///
	title(" ", pos(11))
	
** Generate Figure 4	
*Step 1 - gen child mortality for all africa (minus Cameroon)
harmby aclp_region, by(cowcode)
bysort year: egen cm_af = mean(ChildMortality) if aclp_region==1 & country!="Cameroon"
label var cm_af "Average Child Mortality in Sub-Saharan Africa, excluding Cameroon"
tab cm_af

*gen child mortality for all african oil producers (minus Cameroon)
hist oil_gas_valuepop_2014 if aclp_region==1
hist oilincome if aclp_region==1
tab oilincome if aclp_region==1

*Define oilstate as those where oil and gas production value was non-zero for the years 1980-2010
bysort country: gen af_oil_year = 1 if oilincome>=1 & oilincome!=. & aclp_region==1
tab af_oil_year
list country year if af_oil_year == 1

bysort year: egen cm_af_oil = mean(ChildMortality) if af_oil_year == 1 & country!="Cameroon"
label var cm_af_oil "Average Child Mortality in Sub-Saharan African Oil Producers, excluding Cameroon"
graph twoway (line cm_af_oil year)

***Figure 4 - plots the mortality rate for Cameroon, all other SSA oil producers (see note for definition), and all other SSA
graph twoway (line ChildMortality year if country=="Cameroon" & year>=1980 & year<=2015, xlabel(1980(5)2015, ang(45)) scheme(plotplain)) ///
	(line cm_af year if year>=1980 & year<=2015, ytitle("Child Mortality" " ")) ///  
	(line cm_af_oil year if year>=1980 & year<=2015, ///
	xtitle("") legend(order(1 "Cameroon" 2 "All other SSA" 3 "Other SSA Oil Producers")) ///
	title("Figure 4: Child Mortality Trends in Sub-Saharan Africa, 1980-2015") xsize(8) ///
	note(" " "Oil producers defined as country-years where oil income per capita is at least $1, using data from Ross and Mahdavi (2015)"))	

**re-sort data for following robustness analyses

sort cowcode year	
	
**************************
**************************

*Robustness Analyses - reported in Supplemental Materials File


**1 - examine the linearity of the interaction effect
*first re-estimate main model to split into terciles
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
*moderating variable is 'latent_personalim'
centile latent_personalism if e(sample), centile(33 66)
gen ter = . if e(sample)
replace ter = 0 if latent_personalism<0.3587632 & e(sample)
replace ter = 1 if latent_personalism>=0.3587632 & latent_personalism<0.5460467 & e(sample)
replace ter = 2 if latent_personalism>=0.5460467 & e(sample)

label var ter "Low/Med/High Terciles of Power Concentration"
twoway (sc ln_cm oilincome) (lowess ln_cm oilincome), by(ter, legend(off)) xsize(4.5) scheme(plotplain) xtitle("Log Oil Income Per Capita")

by cowcode: gen lagoil = l.oilincome
by cowcode: gen lagpers = l.latent_personalism
by cowcode: gen laginc = l.inc_comb
by cowcode: gen lagpop = l.log_pop
by cowcode: gen lagreg = l.v2x_regime

interflex ln_cm lagoil lagpers laginc lagpop lagreg ldv if year>=1980, fe(cowcode year) cl(cowcode)

*symmetrical (note: must have Stata SE/MP for this to run; Stata IC can not accommodate sufficient right hand side vars)

*first re-estimate main model to split into terciles
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
*moderating variable is 'latent_personalim'
centile oilincome if e(sample), centile(33 66)
tab oilincome if e(sample)
tab oil_gas_valuepop_2014 if e(sample)
gen oilter = . if e(sample)
replace oilter = 0 if oil_gas_valuepop_2014<1 & e(sample)
replace oilter = 1 if oil_gas_valuepop_2014>=1 & oil_gas_valuepop_2014<91 & e(sample)
replace oilter = 2 if oil_gas_valuepop_2014>=91 & e(sample)
twoway (sc ln_cm latent_personalism) (lowess ln_cm latent_personalism), by(oilter, legend(off)) xsize(4.5) scheme(plotplain) xtitle("Autocratic Personalization")
*interflex ln_cm lagpers lagoil laginc lagpop lagreg ldv if year>=1980, fe(cowcode year) cl(cowcode) type(linear)
*Note: above code is commented out incase users do not have Stata MP, which would cause code to hit an error and stop

** 2 - Robustness analysis

 
*Estimate robustness to replacing ROW ordinal regime measure with interval level polyarchy score
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_polyarchy i.year ldv if year>=1980, fe cluster(cowcode)
estimates store Polyarchy
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90%CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(3.5) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11)) ///
	note("Regime type measured as VDem polyarchy") 

*Estimate robustness to controlling for party regime (as a proxy of whether or not leaders need to personalize; maybe they dont have to personalize because there is a party
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.gwf_party l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
estimates store Party
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90% CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(3.5) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11)) ///
	note("Control for Party Autocracies") 
	
*Estimate robustness to excluding Monarchies	
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980 & gwf_monarchy!=1, fe cluster(cowcode)
estimates store NoMonarchy
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90% CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(3.5) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11)) ///
	note("Drop Monarchies") 
	
*Estimate robustness to controlling for "rigorous and impartial public administration" and "political corruption index" and "conflict" (first two from VDem, conflict from QoG/UCDP/PRIO)
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime l.v2clrspct l.v2x_corr l.conflict i.year ldv if year>=1980, fe cluster(cowcode)
estimates store CapCorrConflict
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90% CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(3.5) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11)) ///
	note("Control for Corruption, Administrative Capacity, & Conflict") 

***Group 1 - plot of coefficients
*Plot these regression coefficients
coefplot Polyarchy Party NoMonarchy CapCorrConflict, drop(_cons *.year ldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income" L.gwf_party = "Party Autocracy" ///
	L.v2x_polyarchy democracy index = "Polyarchy Index" L.v2clrspct = "Administrative Capacity" ///
	L.v2x_corr = "Corruption" L.conflict = "Internal Conflict") ///
	title("Fixed Effect Regression Coefficients and 95% and 90% CI's") ///
	xsize(8) ///
	note(" " "DV is (log) Child Mortality. Sample includes all autocratic regime years 1980-2010." ///
		"All IV's lagged one year. Lagged DV, constant, year, and country fixed effects included" ///
		"but not reported. N/clusters = 2,135/103 for models Party and Capacity/Corruption/Conflict;" ///
		"2,021/99 for model Polyarchy; 1,903/103 for model NoMonarchy. Standard errors clustered on country." ///
		"See text for full sample, variable, and model description." ///
		)
		
		
** 3 - Robustness to using a differenced DV model	
*without the LDV
xtreg d.ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year if year>=1980, fe cluster(cowcode)
estimates store NoLDV
*with LDV (must generate a different lag dv measure, since the dv is now a differenced measure)
by cowcode: gen difflagdv = l.d.ln_cm
xtreg d.ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year difflagdv if year>=1980, fe cluster(cowcode)
estimates store LDV
coefplot NoLDV LDV, drop(_cons *.year LD.ln_cm difflagdv) xline(0) xsize(8) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income") ///
	title("Fixed Effect Regression Coefficients and 95% and 90% CI's" "1 year Differenced DV") ///
	note(" " "DV is annual change in (log) Child Mortality. Sample includes all autocratic" ///
		"regime years 1980-2010. All IV's lagged one year." ///
		"Constant, year, and country fixed effects included but not reported." ///
		"N = 2,135, Clusters = 103; Standard errors clustered on country." ///
		"See text for full sample, variable, and model description." ///
		)

*Models are basically the same, so report ME plots for the no LDV model only, for now
xtreg d.ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year if year>=1980, fe cluster(cowcode)
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90% CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	xsize(8) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	title("Marginal Effect of Oil Income on Child Mortality", pos(11))
	

** 4 - Estimate robustness to RE models with fixed features

xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv l.aclp_region l.ht_colonial l.al_ethnic if year>=1980, cluster(cowcode)
estimates store RandomEffects
*Calculate and plot marginal effects		
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90% CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	xsize(8) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of power concentration, in percent", axis(2)) ///
	legend(order(2 "Marginal Effect" 3 "Distribution of Autocratic Personalization")) ///
	yline(0)) ///
	title("Marginal Effect of Oil Income on Child Mortality") ///
	note("Random Effects model, controlling for region, colonial heritage, and ELF") legend(off)

coefplot RandomEffects, drop(_cons *.year ldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income" L.gwf_party = "Party Autocracy" ///
	L.v2x_polyarchy democracy index = "Polyarchy Index" L.v2clrspct = "Administrative Capacity" ///
	L.v2x_corr = "Corruption" L.aclp_region = "Region" L.ht_colonial = "Colonial Heritage" L.al_ethnic = "ELF") ///
	title("Random Effects Regression Coefficients and 95% and 90% CI's") ///
	xsize(8) ///
	note(" " "DV is (log) Child Mortality. Sample includes all autocratic regime years 1980-2010." ///
		"Model replaces country-specific fixed effects with random effects and a measure of region," ///
		"colonial heritage, and ELF. Lagged DV, constant, and year effects included in each model but not reported." ///
		"N/clusters = 2,024/101. Standard errors clustered on country. See text for full sample," ///
		"variable, and description.")			


** 4 - Alternative: estimate the FE model while restricting to panels with the largest 50% of (within-panel) variation in latent_personalism measure)
*run baseline model
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
*generate a panel specific measure of standard deviation of latent_personalism (the sluggish variable)
bysort cowcode: egen pers_std = sd(latent_personalism) if e(sample)
*define the median of these values; above the median indicates top 50% of distribution in terms of within panel variability in sluggish variable latent_personalism
centile pers_std, centile(50)
*re-run baseline FE model limiting to these observations
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980 & pers_std>0.1026103, fe cluster(cowcode)
estimates store FE_limit
*Plot these regression coefficients as Figure 2
coefplot FE_limit, drop(_cons *.year ldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income") ///
	title("Fixed Effect Regression Coefficients and 95% and 90% CI's" " ", pos(11)) ///
	xsize(8) note(" " "DV is (log) Child Mortality. Sample is autocratic regimes (1980-2010) whose within-panel variability" ///
	"in terms of autocratic personalization is in the top half of the sample distribution. All IV's lagged one year." ///
	"Lagged DV, constant, year, and country fixed effects included but not reported. N = 1,066, Clusters = 50." ///
	"Standard errors clustered on country. See text for full sample, variable, and model description." ///
		)

*Calculate and plot marginal effects
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(8) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11))
	

** 5 - No Lags on IVs
xtreg ln_cm latent_personalism oilincome c.latent_personalism#c.oilincome inc_comb log_pop v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
estimates store NoLags
margins, dydx(oilincome) cont at(latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90% CI)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(white) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	xsize(8) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(order(2 "Marginal Effect" 3 "Distribution of power concentration")) ///
	yline(0)) ///
	title("Marginal Effect of Oil Income on Child Mortality") ///
	note("Fixed Effects Model, no lag on IVs") legend(off)

coefplot NoLags, drop(_cons *.year ldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(latent_personalism = "Autocratic Personalization" oilincome = "Log oil income per capita" ///
	inc_comb = "Log income per capita" log_pop = "Log Population" v2x_regime = "Regime Type - VDem" ///
	c.latent_personalism#c.oilincome = "Autocratic Personalization (x) Oil Income" gwf_party = "Party Autocracy" ///
	v2x_polyarchy democracy index = "Polyarchy Index" v2clrspct = "Administrative Capacity" ///
	v2x_corr = "Corruption" ) ///
	title("Fixed Effects Regression Coefficients and 95% and 90% CI's") ///
	subtitle("No Lag on IVs") ///
	xsize(8) ///
	note(" " "DV is (log) Child Mortality. Sample includes all autocratic regime years 1980-2010." ///
		"Model drops the one-year lag on all independent variables reported in main text." ///
		"Lagged DV, constant, and year effects included in each model but not reported." ///
		"N/clusters = 2,050/101. Standard errors clustered on country. See text for full" ///
		"sample, variable, and description." ///
		)	
			
** 6 - Robustness to Autocratic Instittuions Critique
gen aleg = .
replace aleg = 1 if lparty==2
replace aleg = 0 if lparty<2
tab aleg, mis
label var aleg "Autocratic Legislature"

pwcorr aleg latent_personalism if latent_personalism!=., sig obs
tab aleg if latent_personalism!=.
/*shows that amongst autocracies with parties, about half hanve a legislature
But amongst those with legislatures, 
*/

**full original model, along with the interaction with legislature 
xtreg ln_cm l.aleg l.oilincome c.l.aleg#c.l.oilincome l.latent_personalism c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980, fe cluster(cowcode)
gen insamp=1 if e(sample)
estimates store leg

margins, dydx(l.oilincome) cont at(l.aleg = (0(1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Legislatures") ///
	xlab(-.25 " " 0 "No Legislature" .2 " " .4 " " .6 " " .8 " " 1 "Elected Legislature" 1.25 " ", notick ) ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist aleg if e(sample) , ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	barwidth(.25) ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%", axis(2)) ///
	ytitle("Distribution of Autocratic Legislatures, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(4) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11))	

margins, dydx(l.aleg) cont at(l.oilincome = (0(1)11)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Oil income per capita (log)") ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist oilincome if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50" , axis(2)) ///
	ytitle("Distribution of oil income, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(4) ///
	title("Marginal Effect of Autocratic Legislature on Child Mortality" " ", pos(11))

**Repeat the process, but dropping the power concentration / personalis measure and interaction
*Still need to limit to autocracies (see if 'insamp==1' at end of code, corresponding to above model)
xtreg ln_cm l.aleg l.oilincome c.l.aleg#c.l.oilincome  l.inc_comb l.log_pop l.v2x_regime i.year ldv if year>=1980 & insamp==1, fe cluster(cowcode)
estimates store leg_noper

*plot both models on single coef plot
coefplot leg leg_noper, drop(_cons *.year ldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.aleg = "Elected Legislature" cL.aleg#cL.oilincome = "Elected Legislature (x) Oil Income" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income") ///
	title("Fixed Effect Regression Coefficients and 95% and 90% CI's" " ", pos(11)) ///
	legend(order( 3 "Includes Autocratic" "Personalization"  6 "Excludes Autocratic" "Personalization")) ///
	xsize(8) note(" " "DV is (log) Child Mortality. Sample includes all autocratic regime years 1980-2009. All IV's lagged one year." ///
	"Lagged DV, constant, year, and country fixed effects included but not reported. N = 2,024, Clusters = 101." ///
	"Standard errors clustered on country. See text for full sample, variable, and model description." ///
	 ) 
	
margins, dydx(l.oilincome) cont at(l.aleg = (0(1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Legislatures") ///
	xlab(-.25 " " 0 "No Legislature" .2 " " .4 " " .6 " " .8 " " 1 "Elected Legislature" 1.25 " ", notick ) ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist aleg if e(sample) , ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	barwidth(.25) ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%", axis(2)) ///
	ytitle("Distribution of Autocratic Legislatures, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(4) ///
	title("Marginal Effect of Oil Income on Child Mortality" " ", pos(11))	

margins, dydx(l.aleg) cont at(l.oilincome = (0(1)11)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Oil income per capita (log)") ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist oilincome if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%", axis(2)) ///
	ytitle("Distribution of oil income, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	xsize(4) ///
	title("Marginal Effect of Autocratic Legislature on Child Mortality" " ", pos(11))


** 7 - Robustness to clustering on 'gwf_caseid' (GWF autocratic regime id)
*must run the below line or else there will be repeated time values within panels
drop if latent_personalism == .
tsset gwf_caseid year
*generate new lag DV measure, sorted by gwf_caseid
by gwf_caseid: gen newldv = l.ln_cm
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year newldv if year>=1980, fe cluster(gwf_caseid)
estimates store Standard

*with party dummy
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.gwf_party l.inc_comb l.log_pop l.v2x_regime i.year newldv if year>=1980, fe cluster(gwf_caseid)
estimates store Party

*drop monarchies
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year newldv if year>=1980 & gwf_monarchy!=1, fe cluster(gwf_caseid)
estimates store NoMonarchies

*with differenced dv
by gwf_caseid: gen newdiffldv = l.d.ln_cm
xtreg d.ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year if year>=1980, fe cluster(gwf_caseid)
estimates store DiffDV

*this coefplot call puts the estimates from all four models on the same plot
coefplot Standard Party NoMonarchies DiffDV, drop(_cons *.year newldv newdiffldv) xline(0) scheme(plotplain) levels(95 90) ///
	coeflabels(L.latent_personalism = "Autocratic Personalization" L.oilincome = "Log oil income per capita" ///
	L.inc_comb = "Log income per capita" L.log_pop = "Log Population" L.v2x_regime = "Regime Type - VDem" ///
	cL.latent_personalism#cL.oilincome = "Autocratic Personalization (x) Oil Income" oL.gwf_party = "Party Autocracy") ///
	title("Fixed Effect Regression Coefficients and 95% and 90% CI's" "Cluster SE's by Autocratic Regime") ///
	xsize(8) note(" " "DV is (log) Child Mortality. Sample includes all autocratic" ///
		"regime years 1980-2010, with and without Monarchies. All IV's lagged" ///
		"one year. Constant, year, and country fixed effects included" ///
		"but not reported. Lagged DV included in models 'Full' 'Party'" ///
		"and 'NoMonarchies.' N = 1,970, Clusters = 144 (1,738/135 when" ///
		"dropping monarchies). Standard errors clustered on autocratic regime." ///
		"See text for full sample, variable, and model description." ///
		)
			
*use 'model 1' as the basis for the calculation of ME's (repeat the model in the next line of code)
xtreg ln_cm l.latent_personalism l.oilincome c.l.latent_personalism#c.l.oilincome l.inc_comb l.log_pop l.v2x_regime i.year newldv if year>=1980, fe cluster(gwf_caseid)
margins, dydx(l.oilincome) cont at(l.latent_personalism = (0(.1)1)) atmeans	
marginsplot, level(90) scheme(plotplain) recastci(rspike) xtitle("Autocratic Personalization") ///
	ytitle("Marginal Effect (90%)") ///
	addplot(hist latent_personalism if e(sample), ///
	blcolor(gs10) fcolor(none) ///
	bfcolor(none) ///
	percent ///
	yaxis(2) ///
	xsize(8) ///
	yscale(alt lcolor() axis(2)) ///
	ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%", axis(2)) ///
	ytitle("Distribution of Autocratic Personalization, in percent", axis(2)) ///
	legend(off) ///
	yline(0)) ///
	title("Marginal Effect of Oil Income on Child Mortality", pos(11))

	
**************************
**************************
*generate some descriptive statistics, not just for main sample, but for the universe of possible country-years, and report share of missing data

*define universe as country-years defined as authoritarian by GWF 2018
sum latent_personalism if year>=1980
gen samp_total = 1 if latent_personalism!=. & year>=1980
sum samp_total
tab samp_total, mis
tab year if samp_total==1

*now summarize all variables used in models if samp_total==1, and report N, means, se's, and share missing (N/2290)
sum latent_personalism if samp_total==1
sum ln_cm d.ln_cm latent_personalism oilincome inc_comb log_pop v2x_regime year ///
	v2x_polyarchy gwf_party gwf_monarchy v2clrspct v2x_corr conflict ///
	ht_colonial al_ethnic aleg if samp_total==1
	
	
